Packing Hyperspheres in High-Dimensional Euclidean Spaces 
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Abstract 

We present the first study of disordered jammed hard-sphere packings in four-, five- and six- 
dimensional Euchdean spaces. Using a colhsion-driven packing generation algorithm, we obtain 
the first estimates for the packing fractions of the maximally random jammed (MRJ) states for 
space dimensions d = 4, 5 and 6 to be 0a/_rj — 0.46, 0.31 and 0.20, respectively. To a good 
approximation, the MRJ density obeys the scaling form 4>mrj = ci/2°'+(c2(i)/2'^, where ci = —2.72 
and C2 = 2.56, which appears to be consistent with high-dimensional asymptotic limit, albeit with 
different coefficients. Calculations of the pair correlation function g2{r) and structure factor ^(A;) 
for these states show that short-range ordering appreciably decreases with increasing dimension, 
consistent with a recently proposed "decorrelation principle," which, among othe things, states 
that unconstrained correlations diminish as the dimension increases and vanish entirely in the 
limit d — > oo. As in three dimensions (where (pMRj — 0.64), the packings show no signs of 
crystallization, are isostatic, and have a power-law divergence in g2{r) at contact with power-law 
exponent ~ 0.4. Across dimensions, the cumulative number of neighbors equals the kissing number 
of the conjectured densest packing close to where g2i''') has its first minimum. Additionally, we 
obtain estimates for the freezing and melting packing fractions for the equilibrium hard-sphere 
fluid-solid transition, (pp ~ 0.32 and 0a/ ~ 0.39, respectively, for d = 4, and cpp ~ 0.19 and 
(j)M — 0.24, respectively, for d = 5. Although our results indicate the stable phase at high density 
is a crystalline solid, nucleation appears to be strongly suppressed with increasing dimension. 



* Electronic address: |torquato@electroii.princetoii.e( 
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I. INTRODUCTION 



Hard-sphere systems are model systems for understanding the equihbrium and dynamical 
properties of a variety of materials, including simple fluids, colloids, glasses, and granular 
media. The hard-sphere potential is purely repulsive; it is infinite when two spheres overlap, 
but otherwise zero. Despite the simplicity of the potential, hard-sphere systems exhibit rich 
behavior: they undergo a fluid-solid phase transition and can exhibit glassy behavior. Of 
particular recent interest are (nonequilibrium) disordered jammed packings of hard spheres 
and their statistical and mechanical properties, such as the maximally random jammed 
(MRJ) state , pair correlations j3|, isostaticity [3|, and density fluctuations j^. Such 
ackings have been intensely studied computationally in two and three dimensions 

m 

and in this paper we extend these studies to four, five and six 
dimensions. 

A hard-sphere packing in rf- dimensional Euclidean space M°' is an arrangement of congru- 
ent spheres, no two of which overlap. As in a variety of interacting many-body systems (l2j |. 
we expect studies of hard-sphere packings in high dimensions to yield great insight into the 
corresponding phenomena in lower dimensions. Analytical investi gati ons of hard-spheres can 
be readily extended into arbitrary spatial dimension [isl ll^. Iisl lid 17. 
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23, 



271 l28l . |29| and high dimensions can therefore be used as a stringent testing ground 



for such theories. Along these lines and of particular interest to this paper, predictions have 
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in dis- 



been made about long-wavelength density fluctuations [2^ and decorrelation 
ordered hard-sphere packings in high dimensions. Additionally, the optimal packing of hard 
spheres in high dimensions is also of interest in error-correcting codes in communications 
theory 

Our focus in this paper will be the study of hard-sphere packings in four, five and six 
dimensions. Specifically, we consider both equilibrium packings for d = 4 and d = 5 and 
nonequilibrium packings representative of the maximally random jammed state for d = 4, 
d = 5 and d = 6. 

Equilibrium thermodynamic properties of hard-sphere packings for d = 4 and d = 5 
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3l|. For 



have been studied both computationally and with approximate theories 
the low-density fluid, lower-order virial coefficients, B2, -B3, and B4, are known exactly for 
arbitrary dimensionality Q|. Higher-order virial coefficients have been calculated 
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by Monte Carlo simulation, B ^, B e, and By for both d = A and d = h jl6| and B^ for 
(i = 4 [3], and analytically The pair correlation function for equilibrium fluids has 

been studied and a decrease in ordering with increasing dimension was readily apparent js^ . 
Hard-sphere systems have been shown to undergo a (first-order) fiuid-solid phase transition 
by numerical simulations for 3 < c/ < 5 33| and with approximate theories for d as high as 
50 221 . The freezing points for d = A and d = 'b were estimated numerically to occur at 
packing fractions (pp ~ O.50max = 0.31 and (pp ~ O.40max = 0.19, respectively, and it was 
conjectured that freezing occurs at lower packing fractions relative to close packing as the 
dimension increases j3| • The packing fraction (j) is the fraction of space Mf' covered by the 
spheres, i.e., 

= P^i(i?), (1) 

where p is the number density. 



is the volume of a (i-dimensional sphere of radius i?, and r(a;) is the gamma function 

At sufficiently large densities, the packing of spheres with the highest jamming density 
has the greatest entropy because the free-volume entropy dominates over the degeneracy 
entropy. Therefore, the high-density equilibrium phase corresponds to the optimal packing, 
i.e., maximal density. The densest packing for c? = 3 was recently proven by Hales j3^ to be 
attained by the FCC lattice with packing fraction (pmax = tt/vTS = 0.7404 . . .. The kissing 
number Z, the number of spheres in contact with any given sphere, for the FCC lattice 
corresponds to the maximal kissing number Zmax = 12 for d = 3. One of the generalizations 
of the FCC lattice to higher dimensions is the Dd checkerboard lattice, defined by taking a 
cubic lattice and placing spheres on every site at which the sum of the lattice indices is even 
{i.e., every other site). The densest packing for ci = 4 is conjectured to be the lattice. 
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with packing fraction (pmax = vr^/16 = 0.6168 . . . and kissing number Z = Z^ax = 24 
which is also the maximal kissing number in d = 4 For d = 5, the densest packing is 
conjectured to be the lattice, with packing fraction (f)max = 27r /(30V2) = 0.4652 ... and 
kissing number Z = AO 3u\. For d = 6, the densest packing is conjectured to be the "root" 
lattice Eq, with density (pmax = 37r^/(144v^) = 0.3729 . . . and kissing number Z = 72 
The maximal kissing numbers Z^ax for d = 5 and d = 6 are not known, but have the 
following bounds: 40 < Z^ax < 46 for ci = 5 and 72 < Z^ax < 82 for ci = 6 [s^. In very 



high dimensions, it has been suggested that random packings of spheres might have a higher 
density than ordered packings, enabhng the intriguing possibihty of disordered ground states 
and hence thermodynamic glass transitions see also Ref. 

Equilibrium hard-sphere systems for d = 2 and d = 3 crystallize into ordered packings 
upon densification. However, for (i = 3, it has been found both experimentally Q] and 
computationally [llls, ^ that if the system is densified sufficiently rapidly, the system can be 
kept out-of-equilibrium and can jam in a disordered state. A jammed packing is one in which 
the particle positions are fixed by the impenetrability constraints and boundary conditions, 
despite thermal or mechanical agitation of the particles or imposed boundary deformations or 
loads. Depending on the boundary conditions, different iam ming categories can be precisely 



37 



38, 



39| . The density of disordered 



defined, including local, collective and strict jamming 
collectively jammed hard-sphere packings for = 3 is around (p ~ 0.64 for a variety of 
jacking-generation protocols and has traditionally been called random close packing (RCP) 
3| • However, Ref. jlj| showed that RCP is ill-defined because "random" and "close packed" 
are at odds with one another and the precise proportion of each of these competing effects 
is arbitrary. Therefore, Ref. P] introduced the concept of the maximally random jammed 
(MRJ) state to be the most disordered jammed packing in the given jamming category. This 
definition presupposes an order metric ip can be defined such that ip = 1 corresponds to the 
most ordered [i.e., crystal) packing and = corresponds to the most disordered packing, 



in which there are no spatial correlations. Figure Q from Ref. [Ij shows where MRJ lies on 
a schematic diagram of the space of jammed packings in the density-disorder (p-ip plane. 
In this paper, we numerically study MRJ packings of hard spheres for d = 4, 5 and 6 



that are at least co 
of the MRJ states 



lectively jammed and report the first estimates of the packing fractions 
ll in these dimensions to be (pMRj — 0.46, 0.31 and 0.20, respectively. 
We find that short-range ordering exhibited by g2{r) and S{k) appreciably diminishes with 
increasing dimension, consistent with a recently proposed "decorrelation principle" stating 
that unconstrained spatial correlations vanish asymptotically in high dimensions and that 
the n-particle correlation function Qn for any n > 3 can be inferred entirely from a knowledge 
of the number density p and the pair correlation function (72(1") [3) Ql- ^^^'^ explore 
equilibrium properties, in particular the fiuid-solid phase transition, for d = 4 and d = 5, 
and find a decreased tendency to crystallize with increasing dimension. 

This paper is organized as follows: Section II explains the simulation procedure. Section 



5 




FIG. 1: A highly schematic plot of the subspace in the density-disorder (p — ip plane, where strictly 
jammed three-dimensional packings exist, as adapted from Ref. Point A corresponds to the 
lowest-density jammed packing, and it is intuitive to expect that a certain ordering will be needed 
to produce low-density jammed packings. Point B corresponds to the most dense jammed packing, 
which is also expected to be the most ordered. Point MRJ represents the maximally random 
jammed state. The jamming region in the (p-ip plane will of course depend on the jamming category. 
The gray region is devoid of hard-sphere configurations. 



Ill gives equilibrium results for (i = 4 and (i = 5, Section IV gives results for disordered 
jammed packings for d = 4, 5 and 6, and Section V summarizes and discusses our results. 



II. SIMULATION PROCEDURE 



We use event-driven molecular dynamics and a modified Lubachevsky-Stillinger (LS) al- 



41[, as in Ref. 



42| |. to produce collectively-jammed hard-sphere packings. As in 



gorithm 

Ref. j^^, our algorithm uses periodic boundary conditions applied to a hypercubic cell, 
in which a fundamental cell containing spheres is periodically replicated to fill all of 
Euclidean space. We also use the cell method, in which the computational domain is di- 
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vided into cubic cells and only neighboring cells are checked when predicting collisions for 
a given sphere. Since the number of neighboring cells, as well as the number of spheres 
per cell, increases considerably with increasing dimension, working in high dimensions is 
computationally intensive. Additionally, eliminating excessive boundary effects requires on 
the order of ten sphere diameters per simulation box length, i.e., on the order of = 10*^ 
spheres. Due to the increasing computational load with increasing dimension, we cannot 
yet study d > Q. Implementing the near-neighbor list (NNL) techniques from Ref. 42j |. 
as well as parallelization, are necessary in order to study higher dimensions. Dimension- 
independent C++ codes used to generate the data in this paper can be downloaded at 
http : / / cherrypit . princeton . edu/Packing/C++/ 



Starting from a Poisson distribution of points, the points grow into nonoverlapping 
spheres of diameter D at an expansion rate 7 = dD/dt, while the positions of the spheres 
evolve in time according to Newtonian mechanics, augmented with energy non- conserving 
collisions. Spheres receive an extra energy boost after the collision due to the positive expan- 
sion rate. In practice, the starting configurations for our packing algorithm are low density 
random-sequential-addition packings of spheres [2^. As the density increases, statistics, 
such as pressure, are collected. In the limit 7 — > 0, the system is in equilibrium; for small 
but nonzero 7, the system is in quasi-equilibrium; and for large 7, the system is out of equi- 
librium. Eventually, a jammed state with diverging collision rate is reached. For studies of 
amorphous jammed packings, the expansion must be initially fast to suppress crystallization 
and maximize disorder, but at sufficiently high pressure, the expansion rate must be slow 
enough to allow local particle rearrangements necessary to achieve jamming j^. 



III. EQUILIBRIUM AND METASTABLE PROPERTIES 

The temperature in equilibrium systems of hard spheres is a trivial variable; i.e., it does 
not affect equilibrium configurational correlations, leaving only one independent thermody- 
namic state variable, which can be taken to be either the reduced pressure p = PV/NksT or 
the density 0, related through the equation of state (EOS). Hard-sphere systems undergo a 
(first-order) fiuid-solid phase transition, characterized by a melting point, i.e., the density at 
which the crystal thermodynamically begins to melt, and a freezing point, i.e., the density 
at which the fiuid thermodynamically begins to freeze. Equilibrium properties, such as the 
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melting and freezing points, are studied here using small expansion rates (7 = 10~^ — 10~^) 
and periodic rescaling of the average sphere velocity to one, such that the total change 
in kinetic energy of the system, due to the collisions between growing spheres, was kept 
small. Strictly speaking, a positive growth rate yields nonequilibrium packings but equilib- 
rium packings result as the growth rate tends to zero. The packings were "equilibrated" 
by verifying that orders of magnitude of change in the expansion rate did not change the 
resulting equation of state. In this section we only consider four and five dimensions due to 
(presently) prohibitive computational costs for higher dimensions. 

Figure El shows the reduced pressure p as a function of density for (a) simulations of 
d = 4 systems of spheres placed in a D4 lattice with negative expansion rate 7 = —10^^ and 
(b) simulations of = 5 systems of spheres placed in a lattice with negative expansion 
rate 7 = — 10~^. The pressure initially follows the (lower) crystal branch, until the system 
becomes mechanically unstable and jumps onto the (higher) fluid branch. Also plotted 
is the theoretical prediction of Luban and Michels (LM) for the equation of state jlfil ]. 
which agrees well with our numerical results for the fluid branch for d = 4, but less so 
for d = 5. It is a computational observation that crystals become mechanically unstable, 
giving rise to a sudden jump in pressure, at a density close to the freezing point 

HQ- 
Such "superheating" (undercompression) is most likely due to the difficulty of achieving 

coexistence in finite systems, although we are not aware of a theoretical analysis. From the 

results in Fig. |2l we estimate the freezing points for d = A and d = 5 to he (pp — 0.31 — 0.32 

and (pF — 0.19 — 0.20, respectively. 

The melting points for = 4 and d = 5 can also be estimated from the data in Fig. |21 
Since throughout the coexistence region the fiuid and solid have the same absolute pressure 
P, the melting density can be estimated as the density on the crystal branch with the same 
absolute pressure P as that at the freezing point. The coexistence region is plotted in Fig. |21 
and the melting packing fractions for d = 4 and d = 5 are estimated to be 0a/ — 0.38 — 0.40 
and 0A/ ^ 0.24 — 0.26, respectively. We also observe that the reduced pressure at the freezing 
point is pf — 12 in both d = 4 and d = 5, which agrees with the reduced pressure at the 
freezing point for d = 3, pf — 12.3, obtained from free energy calculations j45|. 

The melting point was also estimated for c? = 4 (higher dimensions are presently too 
computationally demanding) by slowly densifying a system of spheres, initially a fiuid, and 
looking for the onset of partial crystallization, again by monitoring the reduced pressure p as 
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FIG. 2: (Color online) Reduced pressure p as a function of density (p, for a range of system sizes 
(see legend), for (a) d = 4 systems of spheres, initially in a D4 lattice, and negative expansion rate 
7 = — 10~^ and (b) d = 5 systems of spheres, initially in a lattice, and negative expansion rate 
7 = —10^^. N was chosen to make a perfect lattice with periodic boundary conditions, i.e., 
N = {2n)'^/2 for n e Z. Also plotted is the theoretical prediction of Luban and Michels (LM) for 
the equation of state [15i] . Curves for larger system sizes lie farther to the right. 
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a function of density 0. Due to the difficulty of observing coexistence in finite systems and 
the relatively high activation barrier, simulated hard-sphere systems become "supercooled" 
(overcompressed) and nucleation does not occur until the melting density is surpassed. Con- 
sequently, the density at which partial crystallization appears for sufficiently slow expansion 
provides a reasonable estimate for the melting density. Near jamming the reduced pressure 
is asymptotically given by the free- volume equation of state j4q . 



which can be inverted to give an estimate 0j of the jamming density. 

Since the pressure increases very rapidly near jamming, it is more convenient to plot the 
estimated jamming density 4'j{4>) instead of the pressure p{(f)), as shown in Fig.Elfor a system 
of 648 spheres in d = 4. In such a plot, the onset of partial crystallization causes a dramatic 
jump in (j)j{(f)), as the jamming density of the crystal is much higher than the jamming 
density of a disordered packing. The intersection of the curves with the line (f)j{(f)) = (f) gives 
the final jamming density. Sufficiently fast expansion suppresses crystallization and leads 
to packing fractions around 0.45 — 0.47. Slower expansion allows for partial crystallization, 
typically around (pM — 0.38 — 0.39, which is our rough estimate of the melting point, in 
agreement with our estimate from the results in Fig. |21 More accurate estimates can only 
be obtained using free-energy calculations. Since crystallization is a nucleated process, it is 
not surprising that the same expansion rates 7 can crystallize at different packing fractions 
and onto different crystal branches, e.g. 7 = 10"*^ (a) and (b) in Fig. |21 
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FIG. 3: (Color online) Left panel: The estimated jamming packing fraction (j)j as a function of 
density (p for systems of 648 spheres for (i = 4 with various expansion rates (see legend and note 
that there are two samples labeled (a) and (b) for 7 = 10~^). For the curves showing no partial 
crystallization {i.e., 7 = 10~^, 10^^, and 10"''), curves with smaller expansion rates have larger 
peak heights. For the curves that show partial crystallization (i.e. 7 = 10"^ (a and b) and 
10"^), curves with smaller expansion rate lie farther to the left. Right panel: The cumulative 
coordination Z{r) (i.e., the number of contacts) for the perfect D4 lattice and for the partially 
crystallized packings at p > 10^^ obtained for expansion rates 7 = 10~^ and 7 = 10~^. The 
jamming packing fraction for the 7 = 10"^ packing is = 0.511, and the jamming packing fraction 
for the 7 = 10"^ packing agreed up to 12 significant figures with the density of the D4 lattice, 
(f) = 7rVl6 ~ 0.617. 



To determine whether the crystallized packings were forming a D4 lattice, the conjectured 
densest packing in four dimensions, we computed the average cumulative coordination num- 
ber Z{r), which is the average number of sphere centers within a distance r from a given 
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FIG. 4: The estimated jamming packing fraction tpj as a function of density for a system of 
10, 000 spheres for d = 4 with various expansion rates. Curves with smaller expansion rates have 
larger peak heights. The curve labeled "mix" corresponds to the following sequence of expansion 
rates: 7 = 10"^ ^j^^ji ^ ^ ]^o, 7 = lO'^ until p = W^, 7 = 10""^ until p = W^, and 7 = 10"^ until 



sphere center. The inset to Fig. El shows Z{r) for a perfect lattice and for the crystalhzed 
packings with 7 = 10~^ and 7 = 10^^ (corresponding colors represent the same packing). 
The sharp plateaus for the D4 lattice correspond to the coordination shells and the number of 
spheres in the first shell is the kissing number Zmax = 24. The packing shown with 7 = 10~^ 
formed a perfect ^4 lattice. The packing shown with 7 = 10^® partially crystallized with a 
final density of ~ 0.511. 

Figure m shows the estimated jamming packing fraction 0j, as in Fig. El but for a sys- 
tem of 10, 000 spheres, instead of 648 spheres, in four dimensions. In contrast to the 648 
sphere system, there is no sign of partial crystallization for the 10, 000-sphere system. In 
fact, molecular dynamics was performed at packing fractions of ^ 0.38 — 0.42 for 10 mil- 
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lion collisions per sphere and there was no significant drop in pressure indicative of partial 
crystallization. The curves in Figs. El and |3] exhibit a bump around 0g — 0.41, suggesting a 
kinetic transition from the fluid branch to a glassy branch. 
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FIG. 5: (Color online) The estimated jamming packing fraction (pj as a function of packing fraction 
<j) for d = 3. Shown are systems of 4096 spheres with various expansion rates and systems of 10, 976 
spheres placed in an FCC lattice with negative expansion rates 7 = — 10~^, — 10~^, and — 10~^ 
(last three curves). Also plotted are approximations to the equilibrium EOS for the fluid phase, 
the coexistence region, and the crystal phase 47|, as well as the Percus-Yevick (FY) EOS for the 
fluid phase. Compare this figure to the curves shown in Figs. and \^ For the curves showing 
no partial crystallization (i.e., 7 = 32 x 10^^, 64 x 10^^, and 128 x 10~^), curves with smaller 
expansion rates have larger peak heights. For the curves that show partial crystallization (i.e., 
7 = 10~^, 4 X 10~^, and 16 x 10^^), curves with smaller expanion rates lie farther to the left. For 
the melting curves (i.e., 7 = —10"^, —10^^, and —10^^), curves with smaller compression rates lie 
farther to the right. 



Figure El shows the estimated jamming packing fraction 0j for systems of spheres for 
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d = 3 with various positive and negative expansion rates, for comparison with the results 
for d = A and = 5 in Figs. |2l El and |3] The locations of the freezing and melting points 



in d = 3 have been determined from free-energy calculations [45 [and good approximations 
to the EOS for both the fluid and crystal phases are known 43|- Our estimates of the 
freezing and melting points as the densities at the onset of melting of a diluted crystal or of 
partial crystallization of a densified fluid, respectively, compare favorably to the true values 
computed from free-energy calculations in d = 3. The bump around (pc — 0.59, analogous 
to the bump in Fig. HI around (pc — 0.41, is often cited as the approximate location of the 
"kinetic" glass transition j^] • Comparing Figs. 0] and El reveals that the melting point and 
suggested kinetic glass transition are closer for d = 4 than for (i = 3, which is a possible 
reason why there is a lower tendency to crystallize for d = 4 than_for d = 3. Similar results 
have been observed for binary hard disks, a model glass former 



IV. DISORDERED JAMMED PACKINGS 



Packings representative of the maximally random jammed (MRJ) state are produced by a 
combination of expansion rates. The expansion rate must be initially high (compared to the 
average thermal velocity) to suppress crystallization and produce disordered configurations 
that are trapped in the neighborhood of a jammed packing. Near the jamming point, the 
expansion rate must be sufficiently slow to allow for particle readjustments necessary for 
collective jamming. Figure |31 shows the final jamming packing fractions of packings created 
using a variety of expansion rates, as the packing fraction at which the curves intersect the 
line (f)j = (f). We see that by increasing the expansion rate, we attain packings with lower 
jamming packing fractions. 

By comparing Fig. (Hand to the analogous plot for a d = 3 system (Fig. El), where it is 
widely accepted that 0Mi?j — 0.64 — 0.65 |l[ Q]; '^^ estimate the MRJ density for c? = 4 
to be (pMRj — 0.460 ± 0.005. A more accurate calculation of (pMRj demands a better 
theoretical understanding of order metrics and how the expansion rate in the algorithm 
affects the ordering in the produced packings; statistical errors are smaller than the effect of 
the packing-generation protocol. Systematic investigation of different protocol parameters, 
as done for c? = 4 in Fig. IH is currently too computationally intensive in higher dimensions. 
Reasonable estimates of (pMRj for both d = 5 and d = 6 are obtained using the following 
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less computationally intensive procedure. First, the system of spheres is expanded, starting 
from zero initial kinetic energy (T = 0), until it reached a high pressure (say, p = 100 — 
1000). Then the system is slowly expanded (7 = 10"^ — 10~^) and periodically cooled to 
ksT = 1 until a very high pressure (say, p = 10^^) is attained. The resulting packings are 
approximately collectively jammed, as demonstrated by very large relaxation times for the 
pressure during long molecular dynamics runs |Q|. Using this method we estimate the MR J 
density for = 5 to be (pMRj ^ 0.310 ± 0.005 and for = 6 to be (pMRj ^ 0.200 ± 0.01. 

The MRJ packing fractions as well as important equilibrium packing fractions are summa- 
rized in Table HI It is useful to compare the MRJ packings fractions for 3 < c/ < 6 to recent 
estimates of the saturation packing fraction (pg for the random sequential addition (RSA) 
packing of hard spheres obtained by Torquato, Uche and Stillinger 2^ in corresponding 
dimensions, which were shown to be nearly hyperuniform 2^. These authors found that 
(ps = 0.38278 ±0.000046, 0.25454 ±0.000091, 0.16102 ±0.000036 and 0.09394 ± 0.000048 for 
d = 3,4,5 and 6, respectively. The nonequilibrium RSA packing is produced by randomly, 
irreversibly, and sequentially placing nonoverlapping spheres into a volume. As the process 
continues, it becomes more difficult to find available regions into which the spheres can be 
added. Eventually, in the saturation (infinite-time) limit, no further additions are possible, 
and the maximal achievable packing fraction is the saturation value (pg [see Ref. and 
references therein]. As expected, the RSA saturation packing fraction in dimension d is 
substantially smaller than the corresponding MRJ value because, unlike the latter packing, 
the particles cannot rearrange. 

Our estimates for the MRJ packing fraction are compared to a theoretical formula pro- 



posed by Philipse 



2J| for the "random jamming density" (pa, 



0.046^2 + 1.22ci±0.73 
'^'^ ' (5) 



which predicts 03 ~ 0.601, 04 ^ 0.397, 05 ~ 0.249, and 06 ^ 0.152. It is seen that Eq. © 
underestimates MRJ density (puRj in d = 3 and becomes worse with increasing dimension. 
Following Ref. [2^, we obtain a better scaling form by noting that the product 2'^(pMRj for 
3 < c/ < 6 is well approximated by a function linear, rather than quadratic, in d (see Fig. 
in}, i.e., the scaling form for 0a/_rj is given by 

(pMRJ - ^ + (Dj 
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Packing 










fraction 


d = 3 


d = 4 


d = 5 


d = 6 


4'F 


0.494 [23, 451 


0.32 ±0.01* 


0.19 ±0.01* 




(t>M 


0.545 [23, 45j 


0.39 ±0.01* 


0.24 ±0.01* 




(t>MRJ 


0.645 ± 0.005 [2] 


0.46 ± 0.005* 


0.31 ± 0.005* 


0.20 ±0.01* 


4'max 


0.7405... [34] 


0.6169... [30] 


0.4652... [30] 


0.3729... [3^ 



TABLE I: Important packing fractions for d = 3, 4, 5 and 6. These include the equilibrium values 
for the freezing, (pp, melting, (pM, and densest states, (pmax, as well as the nonequihbrium MRJ 
values. The freezing and melting points for d = 6 were not calculated here. *Values computed in 
this work. 




FIG. 6: (Color online) Fit of the data for the product 2'^(j)MR.j to the hnear form © for 3 < d < 6 
with ci = —2.72 and C2 = 2.56. 



where ci = —2.72 and C2 = 2.56. Although the scaling form (jH} applies only in low di- 
mensions such that d > 3, theoretical arguments given by Torquato, Uche and Stillinger 
2^ suggest that the general scaling form © persists in the high-dimensional asymptotic 
limit, albeit with different coefficients ci and C2. In Ref. ||29|, the density lower bound 
4>MRj > {d + 2)/2'^ is derived for MRJ packings in any dimension. This MRJ density lower 
bound yields 0.3125, 0.1875^ 0.109375, 0.0625 for d = 3,4,5 and 6, respectively. We note 
that Parisi and Zamponi [2^ suggest the MRJ density scaling (pMRj ~ {d\ogd)/2'^. 
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A. Pair Correlations 



Our main interest is pair correlations in the jamming limit in four, five and six dimensions. 
We characterize jammed packings statistically using the pair correlation function g2{r) and 
structure factor S{k). The pair correlation function measures the probability of finding 
a sphere center at a given distance from the center of another sphere, normalized by the 
average number density p to go asymptotically to unity at large r; i.e. 



where P(r) is the probability density for finding a sphere center a distance r from an arbitrary 
sphere center, () denotes an ensemble average, and Si(r) is the surface area of a single 
hypersphere of radius r Si(r) = 27r^r^ in d = 4, Si(r) = Svr^r'^/S in d = 5 and 

Si{r) = TT^r^ in d = 6. The structure factor 



is related to the Fourier transform of the total correlation function h{r) = g2{r) — 1. It mea- 
sures spatial correlations at wavenumber k and in particular, large-scale density fluctuations 
at k = 25] . The structure factor can be observed directly via scattering experiments Q| . 

In the jamming limit, the pair correlation function g2{r) consists of a (5-function due to 
sphere contacts and a background part g2{r) due to spheres not in contact: 



where Z is the average kissing number. Figure U\ compares the pair correlation function for 
jammed packings of 10^ spheres in = 3, 4, 5 and 6. Due to periodic boundary conditions, 
g2{r) can only be calculated up to half the length of the simulation box, which limits the 
calculation to r/D ~ 3 for = 6. The well-known split second peak present in ci = 3 
is strongly diminished as the dimension increases, i.e., the amplitude of the split second 
peak decreases and the sharp cusps become rounded with increasing dimension. The split 
third peak present in c? = 3 with considerable structure and two shoulders vanishes almost 
completely in the higher dimensions. The oscillations are strongly damped with increasing 
dimension and the period of oscillations might also decrease slightly with increasing dimen- 
sion; this latter possibility is revealed more vividly in the structure factor through the shift 




(7) 



S{k) = l + ph{k) 



(8) 




(9) 
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r/D 

FIG. 7: The pair correlation function g2{T) for MRJ packings of 10"^ hard spheres for d = 3, 4, 5 
and 6 at the respective densities reported in Table ^ Pair separation is plotted in units of the 
sphere diameter D. [For d = 6, g2{r) was only calculated up to r/D = 3 due to the system size and 
periodic boundary conditions]. The delta- function contribution [c/. Eq. [5] at contact, of course, is 
not shown. The inset shows \h{r)\ = 1^2 (j") — 1| on a logarithmic scale for d = 3, 4 and 5. Each 
curve for g2{r) is obtained from a single packing realization (not time-averaged). Curves for higher 
dimensions are increasingly diminished. 

in the location of the maximum, as we will describe below. The inset to Fig. [7| shows the 
magnitude of the decaying oscillations in h{r) on a semi-log scale. Though at the values 
of r/D shown, up to about half the length of the simulation box, there is still structure in 
addition to the oscillations, especially apparent for d = 3, it appears that the decay rate of 
the oscillations in h{r) does not change significantly with dimension, whereas the amplitude 
of oscillations does. However, further studies with larger r and therefore larger systems are 
needed to obtain more quantitative results. 
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We calculate the structure factor S{k), defined in Eq. |S1 for c? = 4 and d = 5 hj 

S{K) = 1 + 1280 x^h{x)^^^^^dx (10) 







Kx 



and 



S{K) = 1 + 48( 



x'^h{x) 
{Kxf 



sin(i^x) 
Kx 



cos(-ft'x) 



dx. 



(11) 

r/D and 



respectively, where = 7r^pD^/32 for d = A and = n'^pD^ /60 for d = 5, x 
K = kD are the dimensionless radius and wave number, and Ju{x) is the Bessel function of 
order u. We do not calculate the structure factor for = 6 because at present we do not 
have g2{r) over a sufficiently large range of r. 

Following Ref. j^, rather than working directly with g2{x) as in Eq. (jH)), we consider the 
average cumulative coordination Z{x), defined to be the following volume integral of g2{x)'- 



Z{x) — P J Si{x')g2{x')dx' . 



The excess coordination AZ{x), 

AZ{x) 
AZ{x) 



1 + 640 / {x'yh{x')dx' 





1 + 1600 / {x')^h{x')dx', 
Jo 



(12) 

(13) 
(14) 



for d = A and d = 5, respectively, is the average excess number of sphere centers inside a 
spherical window of radius x centered at a sphere, compared to the ideal gas expectations, 
160x^ for d = 4 and 320x^ in c? = 5. We can rewrite Eq. (jHJ in terms of AZ{x) using 
integration by parts to get 



SiK) 



and 



S{K) 



AZ{x] 



AZ{x)^^l^dx 
dx Kx 

siia{Kx) cos{Kx) 



dx 



{Kx) 



(Kx) 



dx. 



(15) 



(16) 



for d = 4 and d = 5, respectively. Note that accurate evaluations of the integrals of AZ{x) 



require extrapolations of its 
damped oscillating function 



^^e-x tail behavior, for which we have used an exponentially- 

FigurelSl shows S{k) for jammed packings of 10^ spheres in three, four and five dimensions. 
Qualitatively, S{k) is somewhat similar for d = 3, 4, and 5. However, with increasing 
dimension, the height of the first peak of S{k) decreases, the location of the first peak 
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FIG. 8: The structure factor S{k) for jammed packings of 10^ spheres for d = 3, 4 and 5 at the 
respective densities reported in Table d Inset: A comparison for d = 4 of S{k) for a jammed 
packing and for a fluid near the freezing point {(f) ~ 0.31). Each curve for S{k) is obtained from a 
single packing realization (not time averaged). 



moves to smaller wavelengths, and the oscillations become damped. The width of the first 
peak also increases with increasing dimension, which could indicate that the correlation 
length decreases with increasing dimension. The inset to Fig. |S1 shows S{k) for a jammed 
packing and a fluid near the freezing point in four dimensions. The relation between the 
structure factor for the fluid and jammed packing is strikingly similar to what is found 
for d = 3, except that the peaks of both curves for c? = 4 appear scaled down relative to 
d = 3. Overall, our results for both g2{r) and S{k) are consistent with a recently proposed 
"decorrelation" principle j^]. We note that simil ar p air decorrelations are observed for RSA 
packings as the dimension increases up to c? = 6 [2^ . 

It is of interest to determine whether infinite-wavelength density fluctuations S{k = 0) 
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vanish; systems with this property are called "hyperuniform" [25|. For equilibrium fluids 
and crystals, S{k = 0) is proportional to the isothermal compressibility and therefore must 
be positive. As for d = 3, S{k) for d = 4 appears to go to zero faster near the origin for the 
jammed packing than for the fluid. However, we cannot reliably determine whether S{k) 
vanishes at the origin because our calculation of S{k) for small k involved an extrapolation 
of the large-x tail of AZ{x). Nevertheless, using larger system sizes of one million spheres, 
saturated 50| MRJ packings for d = 3 have been shown to be hyperuniform to a high 
accuracy j4| and the comparison of d = A and d = 5 to d = 3, shown in Fig. suggests that 
MRJ packings for d = 4 and d = 5 are also hyperuniform. 

B. Isostaticity 

We study the near-contact contribution to g2{r), i.e., interparticle distances r that are 
very close to the sphere diameter D, using the cumulative coordination number Z{x), where 
as before x = r/D is the dimensionless radius and x — 1 is the dimensionless interparticle 
gap. Figure IHl shows Z(x) for jammed packings of 10, 000 spheres for d = 4 and d = 5 
with rattlers removed [51 1. The plateaus at Z = 8 in Fig. |H1 (a) and Z = 10 in Fig. El (b) 
show that both packings are isostatic. Isostatic packings are jammed packings which have 
the minimal number of contacts necessary for collective jamming. For spheres, this occurs 
when the number of degrees of freedom equals the number of contacts (or constraints); 
each d-dimensional sphere has d degrees of freedom, and hence the mean number of contacts 
experienced by a sphere necessary for jamming is 2d, since each contact involves two spheres. 

Packings produced by the LS algorithm almost always contain a nonzero fraction of 
"rattlers", which are spheres trapped in a cage of jammed neighbors, but free to move 
within the cage. We find approximately ~ 1% rattlers for d = 4 and ~ 0.6% rattlers for 
= 5, as compared to ~ 2 — 3% rattlers for = 3 Jsl- Rattlers can be identified as having 
less than the required d + 1 contacts necessary for local jamming and are removed to study 
the jammed backbone of the packing, which we focus on in this section. 

The insets to Fig. 121 (a) and (b) show Z{x) —2d, along with a power-law fit for intermediate 
interparticle gap x — 1, 

Z{x) = 2 + Zo{x - 1)", (17) 
where Z = 2d. Since the packings are generally slightly subisostatic, we apply a small 
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correction (< 0.1%) to the isostatic prediction of 2d by using the midpoint of the apparent 
plateau in Z{x). The best-fit exponent is a ~ 0.6 in both d = A and = 5, in agreement with 
that found for = 3 j^. The coefficients of the power law, ~ 11 in c? = 3, ~ 24 for 
d = 4, and Zq ^ 40 for d = 5 are close to the corresponding kissing numbers of the densest 
packings, Z = 12 for = 3, Z = 24 for = 4, 40 < Z < 46 for c/ = 5 and 72 < Z < 80 
for d = 6. Motivated by this observation, we measured the value of the gap x — 1 at which 
the cumulative coordination Z{x) equals the kissing number of the densest packing to be: 
X — 1 ^ 0.35, 0.34, 0.31 — 0.36 and 0.33 — 0.36 in c? = 3, 4, 5 and 6, respectively, which we can 
define to be the cutoff for the near-neighbor shell. This definition produces results similar 
to that of the more common definition of the cutoff for the near-neighbor shell as the value 
of the gap X — 1 at the first minimum in g2, which occurs at x — 1 ^ 0.35, 0.32, 0.30 and 
0.28 in d = 3, 4, 5 and 6, respectively. It is also interesting to observe that the power-law fit 
to Z{x) is good over a rather wide range of gaps, almost up to the first minimum in g2. We 
should, however, emphasize that the minimum of g2 is not very precisely defined, especially 
due to decorrelation in high dimensions, and the choice of the gap at the minimum of g2, 
or at which Z{x) equals the kissing number of the densest packing, as a special point is 
somewhat arbitrary and not theoretically justified at present. 



V. DISCUSSION 



We have presented the first numerical results characterizing random jammed hard-sphere 
packings in four, five and six dimensions. We find disordered packings, representative of the 
maximally random jammed state, to be isostatic and have packing fractions (pMRj — 0.46, 
(pMRj — 0.31 and 0a/_rj — 0.20 for = 4, 5 and 6, respectively. For equilibrium sphere 
packings, we estimate the freezing and melting packing fractions for the fiuid-solid transition 
in four dimensions to be (pp ~ 0.32 and 0m — 0.39, respectively, and in five dimensions to 
be (pF ^ 0.19 and 0m ~ 0.24, respectively. Additionally, a signature characteristic of the 
kinetic glass transition is observed around (pc — 0.41 for d = 4. We observe a significantly 
lower tendency to crystallize for d = 4 than in d = 3, which is likely due to the closer 
proximity of the melting and kinetic glass transition densities for c? = 4 ^j. 

We find that in high dimensions the split-second peak in the pair correlation function g2, 
present for d = 3, gets dramatically diminished and oscillations in both g2 and the structure 
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FIG. 9: The near-contact cumulative coordination Z{x) [c.f. Eq. for IC^-sphere MRJ packings 
for (i = 4 (a) and for d = 5 (b), with rattlers removed. The inset shows Z{x) on a log- log scale 
along with power-law fits for intermediate interparticle gap x — 1 beyond contact. 10^-sphere MRJ 
packings in d = 5 with final expansion rates of 7 = 10"'^ give similar results; such packings with 
final expansion rates of 7 = 10~^ are (presently_) too computationally expensive. Compare these 
plots to the equivalent results for d = 3 in Ref.2^ [c.f. Fig. 8]. 



factor S{k) get significantly danipened. Tliese findings are consistent with a recently pro- 
posed "decorrelation principle" j27|, stating that unconstrained spatial correlations vanish 
asymptotically in the high- dimensional limit and that the n-particle correlation function 
Qn for any n > 3 can be inferred entirely from a knowledge of the number density p and 
the pair correlation function (72(5"). Accordingly, in this limit the pair correlation function 
g2{r) would be expected to retain the delta-function contribution from nearest-neighbor 
contacts, but the extra structure representing unconstrained spatial correlations beyond a 
single sphere diameter would vanish. Figures [7| and |H1 show dramatically the decorrelation 
principle already taking effect in four, five and six dimensions. We note that decorrelation 
principle is also apparent in the same dimensions for RSA packings Q|. 

One should not be misled to believe that the decorrelation principle is an expected "mean- 
field" behavior. For example, it is well known that in some spin systems correlations vanish 
in the limit d 00 and the system approaches the mean-field behavior. While this idea 
has meaning for spin systems with attractive interactions, hard-core systems, whose total 
potential energy is either zero or infinite, cannot be characterized by a mean field. Mean- 
field theories are limited to equilibrium considerations, and thus do not distinguish between 
"constrained" and "unconstrained" correlations because, unlike us, they are not concerned 
with non-equilibrium packings of which there an infinite number of distinct ensembles. The 
decorrelation principle is a statement about any disordered packing, equilibrium or not. 
For example, contact delta functions are an important attribute of non-equilibrium jammed 
disordered packings and have no analog in equilibrium lattice models of any dimension. 
The decorrelation principle is also justified on the basis of a rigorous upper bound on the 
maximal packing density in high dimensions [27j, which has no counterpart in mean- field 
theories. 

A particularly interesting property of jammed hard-sphere packings is hyperuniformity, 
the complete suppression of infinite wavelength density fluctuations, i.e., the vanishing of 



the structure factor S{k) as — > 0. It has 
strictly-jammed packings are hyperuniform [2] 



)een recently conjectured that all saturated 
and calculations of the structure factor near 



k = for d = 3 using one million particle systems have strongly suggested that MRJ 
packings for d = 3 are indeed hyperuniform ^ . Though the system sizes used in this paper 
were too small to probe such large-scale density fluctuations without relying on dubious 
extrapolations, our numerical results for the structure factor for d = 4 and ci = 5, as shown 
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in Fig. ISl are consistent with hyperuniformity. 

As in three dimensions, disordered jammed sphere packings show no signs of crystaUiza- 
tion, are isostatic, and have a power-law divergence in g2{r) at contact. Interestingly, all 
three dimensions (3, 4 and 5) share the same power law exponent 1 — a ~ 0.4 when rattlers 
are removed, and show the first minimum of g2{r) close to where the cumulative coordina- 
tion Z{r) equals the kissing number of the densest lattice packing. Such a relation between 
the kissing numbers of the densest packings and MRJ packings for d = 3, 4, 5 and 6, if not 
coincidental, is very surprising and may be a consequence of the geometrical structure of 
MRJ packings. It suggests that disordered packings might be deformed crystal packings, in 
which the true contacts are deformed into near contacts, and only the minimal number of 
contacts necessary for jamming is preserved. This interpretation is to be contrasted with 
the usual interpretation of disordered packing in = 3 in terms of tetrahedral or icosahe- 
dral packings, without relation to the crystal (FCC) packing. The former interpretation is 
similar to the one of the MRJ state for binary hard disks as a random partitioning of the 
monodisperse triangular crystal into "small" and "large" disks, i.e., a deformed monodis- 
perse triangular disk crystal in which a randomly chosen fraction of the particles have grown 
in size, as proposed in Ref. Q]. 

It is important to point out that hard-sphere packings behave rather differently in two 
dimensions than in three and higher dimensions. For d = 2, jammed hard-sphere systems 
are polycrystalline and there is a very weak, nearly continuous fluid-solid phase transition. 
Hence, there is no glassy behavior for d = 2 and consequently no amorphous jammed pack- 
ings. Glassy behavior, due to geometrical frustration arising from the inconsistency of local 
optimal packing rules and global packing constraints, first appears in three dimensions . 
It is likely that geometrical frustration generally increases with dimension, consistent with 
our observation that nucleation is suppressed with increasing dimension. 

Computational costs rise dramatically with increasing dimension and theoretical under- 
standing based on observations in moderate dimensions is necessary. We believe that the 
numerical results presented in this work provide tests and motivations for such theories. 
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